//----------------------------------------------------------------------------
// M A R I T I M E  R E S E A R C H  I N S T I T U T E  N E T H E R L A N D S
//----------------------------------------------------------------------------
// Copyright (C) 2008 - MARIN - All rights reserved - http://www.marin.nl
//----------------------------------------------------------------------------
// Program    : mermaid
// Module     : vbm
// File       : RrbElwinSolver.cpp
//----------------------------------------------------------------------------

#ifdef _WINDOWS
#define NOMINMAX
#endif

#include <math.h>
#include <assert.h>
#include "RrbElwinSolver.h"
#include "FloatUtils.h"
#include "defines.h"
#include "Timer.h"

#include <string>
#include <stdio.h>
#include <fstream> // voor printf
#include "assertplus.h"

#ifdef _WINDOWS
#pragma warning(disable: 4244) // conversion from 'double' to 'float', possible loss of data
#endif

#ifdef EWLAPACK
// gebruik lapack ipv eigen routines
#include <lapackpp.h> // Lapack++ library // see NopDefSolver.cpp
#include <sybfd.h> // for LaSymmBandMatFactorizeIP()
#endif

#define COLUMNS  (unsigned)(Nx1 + 3)
#define ROWS     (unsigned)(Nx2 + 3)

// extra parameter voor maximum level toegevoegd (rrb - k)

//  March 2009, Elwin van 't Wout
//      added comments
//      changed the order in the CG-iteration
//      added maxIterations
//      changed termination criterium
//      declaration of variables at begin of routine
//      renamed some variables:
//          alf -> alpha; bet -> beta; telmet -> iteration; pq -> sigma; rhobew -> stopvalue;
//          B -> r; X -> x; Y -> y; Az -> q;
//          m_cne -> p; m_cse -> z; m_cnw -> y
//      changed the numbering of levels niv
//      static void matvec() -> private void matvec()
//      matvec(cc,cs,cw,cne,cse,cnw,Nx1,Nx2) -> matvec(p,q,y)
//      void solver() -> void precsolve()
//      precsolve() -> precsolve(z)
//      m_cne, m_cnw, m_csw, m_csw become local in precsolve()

//  =================================================================== 
//          List of symbols
//  =================================================================== 
//
//  Nx1, Nx2
//      Number of unknowns in x- resp. y-direction.
//
//  tolerance, stopcriterium
//      tolerance: The tolerance of the iterative solver. This is an input variable.
//      stopcriterium: The value of the termination criterium.
//
//  maxIterations
//      Maximum number of CG-iterations allowed.
//
//  maxLevel, m_nivmax
//      maxLevel: Maximum number of RRB-levels as given as input.
//      m_nivmax: Maximum number of RRB-levels as calculated with (<= maxLevel).
//
//  cc,cs,cw
//      The center, south and west stencils of the original matrix.
//
//  ccc,cww,cee,cnn,css
//      The stencils of the RRB-preconditioner.
//
//  cne,cnw,cse,csw
//      Used in the preconditioner for the nine-point stencils.
//
//  x
//      Solution and iterand during CG-iteration.
//
//  r,z
//      r: Right hand side and residual during CG-iteration.
//      z: Preconditioned residual during CG-iteration.
//
//  p,q
//      p: Search vector during CG-iteration.
//      q: Preconditioned search vector during CG-iteration.
//
//  iteration, alpha, beta, rho, rhoold, rhonew, sigma
//      Parameters in CG-iteration
//      iteration: The number of iteration.
//      alpha: rhonew/sigma.
//      beta: rhonew/rhoold.
//      rho: <r,r>, the Euclidean norm of the residual.
//      rhoold, rhonew: <r,z>, the A-norm of the residual with the preconditioner as matrix.
//      sigma: <p,q>
//
//  fac
//      Parameter to raise the diagonal elements. Similar to Gustafsson's modification.
//      Because the matrix is strictly diagonally dominant, this is not necessary.
//
//  st, tw, niv, fac
//      Variables for RRB-iteration.
//      niv: level of RB-orderings.
//      st: starting grid point of RB-ordering
//      tw: the mesh width at a certain level
//
//  m_dimS, m_bandwith
//      m_dimS: Dimension of the last Schur complement.
//      m_bandwith: Bandwith of the last Schur complement.
//
//  m_Chol
//      Array with the elements of the Cholesky matrix. Only the diagonals are stored.
//

const int maxIterations = 1000; ///< maximum number of CG-iterations
//const int maxLevel = 9;

/// constructor

RrbElwinSolver::RrbElwinSolver(Array2D<REAL> *pcc,
                     Array2D<REAL> *pcs,
                     Array2D<REAL> *pcw,
                     int Nx1,
                     int Nx2,
                     REAL tolerance,
                     int maxLevel)
: base(Nx1, Nx2, tolerance)
, m_pcc(pcc)
, m_pcs(pcs)
, m_pcw(pcw)
, m_ccc(COLUMNS, ROWS)
, m_css(COLUMNS, ROWS)
, m_cww(COLUMNS, ROWS)
, m_cee(COLUMNS, ROWS)
, m_cnn(COLUMNS, ROWS)
, m_Chol()
#ifdef EWLAPACK
, m_LaChol()
#endif
, m_dimS()
, m_bandwith()
, m_nivmax(maxLevel)
{
    preconditioner();
}

/// destructor

RrbElwinSolver::~RrbElwinSolver()
{
}

/// calculates preconditioner
/** Deze subroutine maakt een LDL'-ontbinding van een symmetrische, pos. sem. 
    definiete matrix volgens de RRB-methode. \n
    L' en D worden opgeslagen in stencil-vorm. 
    L' wordt in opgeslagen in de arrays cnn, css, cee en cww. \n
    Omdat het nieuwe Schurcomplement telkens gerepresenteerd kan worden als een 
    negen-punts stencil, hebben we tijdens de incomplete ontbinding als kladarrays 
    nodig: cse, csw, cne en cnw. \n
    Het laatste stuk wordt exact gedaan. Het laatste Schur-complement heeft 
    grootte sg en bandbreedte bsm. Een incomplete ontbinding van het laatste 
    Schur-complement wordt opgeslagen in s.  */
/* EW: Nu wordt de preconditioner in een matrix over het hele domain genomen. Maar je hebt alleen
 * de rode punten nodig, omdat de zwarte punten van het eerste niveau niet in het Schur-complement komen.
 * Hierdoor kunnen de arrays ccc,cee,cnn,cww,css,cne,cnw,csw,cse de helft kleiner. Maar dan
 * wel minder leesbare loops. */

/* EW: Het maximaal aantal levels tijdens de RRB-preconditionering was eerst een 'const int',
 * nu is dit een invoerparameter. De vierde parameter van lin_wacu.exe (in de commandoregel) geeft
 * deze parameter aan. */
void RrbElwinSolver::preconditioner()
{

    Array2D<REAL> &cc = *m_pcc; // center stencil original system
    Array2D<REAL> &cs = *m_pcs; // south stencil original system
    Array2D<REAL> &cw = *m_pcw; // west stencil original system

    Array2D<REAL> cne(cc.dim1(), cc.dim2()); // north-east stencil during RRB-procedure
    Array2D<REAL> cnw(cc.dim1(), cc.dim2()); // north-west stencil during RRB-procedure
    Array2D<REAL> csw(cc.dim1(), cc.dim2()); // south-west stencil during RRB-procedure
    Array2D<REAL> cse(cc.dim1(), cc.dim2()); // south-east stencil during RRB-procedure

#ifndef EWLAPACK
    REAL sumc; // dummy variable during the Cholesky decomposition
#endif
    REAL omega; // parameter of the lumping procedure
    omega = 1.0; // default 1.0

    int ii, jj; // counting variables
    int st, tw, niv; // dummy variables during RRB-procedure
    int NSx1, NSx2; // number of nodes in x- and y-direction of last Schur complement

    /* EW: waarom zijn dit allemaal member variabelen en niet lokale. Hierin staat het stencil
     * van de preconditionering, ook nodig tijdens CG-iteratie. */
    // Als eerste stap de matrix A kopieren naar ccc, cww, enz. 
    m_ccc = cc;
    m_css = cs;
    m_cww = cw;

    /* EW: in principe kan je symmetrie gebruiken en zodoende heb je cnn en cee niet nodig,
     * maar dan wordt de code een nog grotere boekhouding. */
    // EW: de loop gaat niet over de hele array !?
    for (jj = 0; jj <= m_Nx2 + 1; ++jj)
    {
        for (ii = 0; ii <= m_Nx1 + 1; ++ii)
        {
            m_cnn[ii][jj] = cs[ii ][jj + 1];
            m_cee[ii][jj] = cw[ii + 1][jj ];
        }
    }

    cne = 0.0;
    cnw = 0.0;
    cse = 0.0;
    csw = 0.0;

    // The number of levels should be at least two.
    // EW: in theory, it can also be one, which yields a complete Cholesky decomposition, however,
    // because the CG will be called anyway, problems occur, probably because a parameter becomes zero.
    if (m_nivmax < 2)
    {
        m_nivmax = 2;
    }

    // Be sure the maximum level has some nodes in it.
    /* EW: waarschijnlijk hoeft dit niet, want anders heeft het maximale niveau
     * alleen maar (0,0) en dan gaat het alsnog goed. */
    //while ( pow(2.0,m_nivmax-1) > std::max(m_Nx1,m_Nx2) )
    while (pow(2.0, m_nivmax - 1 - 1) > std::max(m_Nx1, m_Nx2))
    {
        --m_nivmax;
    }

    st = 1;
    tw = 2;
    niv = 1; // level in the RRB-iteration

    /* Start the RRB preconditioner iteration. Proceed till maximum level nivmax
     * or as long as the grid allows. */
    while (niv < m_nivmax && (st < m_Nx1 + 2 || st < m_Nx2 + 2)) // EW: laatste niet meer nodig
    {
        //                                 cnw cnn cne 
        // Ga uit van een 9-punts stencil  cww ccc cee 
        //                                 csw css cse 
        // Doorloop eerst een schaakbord met daarin cnw, cne, ccc, csw en cse 
        // als zwarte punten in S11. Lump cnw, cne, csw en cse op de rest. 
        // Dit geeft een 5-punts stencil en diagonaalmatrix S11. 
        //                                     cnn 
        //                                 cww ccc cee 
        //                                     css 
        if (niv > 0)
        {
            // Several lumping strategies

            // 1) similar to RIC
            for (jj = st; jj <= m_Nx2 + 2; jj += tw)
            {
                for (ii = 0; ii <= m_Nx1 + 2; ii += tw)
                {
                    m_ccc[ii][jj] = m_ccc[ii][jj] + omega * (
                                                             cne[ii][jj] +
                                                             cnw[ii][jj] +
                                                             cse[ii][jj] +
                                                             csw[ii][jj]);
                }
            }
            for (jj = 0; jj <= m_Nx2 + 2; jj += tw)
            {
                for (ii = st; ii <= m_Nx1 + 2; ii += tw)
                {
                    m_ccc[ii][jj] = m_ccc[ii][jj] + omega * (
                                                             cne[ii][jj] +
                                                             cnw[ii][jj] +
                                                             cse[ii][jj] +
                                                             csw[ii][jj]);
                }
            }
            /*
                        // 2) similar to RIC, but on the outer elements
                        for (jj = st; jj <= m_Nx2 + 2; jj += tw)
                        {
                    for (ii = 0; ii <= m_Nx1 + 2; ii += tw)
                            {
                                m_cee[ii][jj] += omega * (cse[ii][jj] + cne[ii][jj])/2.0;
                                m_cnn[ii][jj] += omega * (cnw[ii][jj] + cne[ii][jj])/2.0;
                                m_cww[ii][jj] += omega * (csw[ii][jj] + cnw[ii][jj])/2.0;
                                m_css[ii][jj] += omega * (csw[ii][jj] + cse[ii][jj])/2.0;
                    }
                    }
                    for (jj = 0; jj <= m_Nx2 + 2; jj += tw)
                        {
                    for (ii = st; ii <= m_Nx1 + 2; ii += tw)
                        {
                                m_cee[ii][jj] += omega * (cse[ii][jj] + cne[ii][jj])/2.0;
                                m_cnn[ii][jj] += omega * (cnw[ii][jj] + cne[ii][jj])/2.0;
                                m_cww[ii][jj] += omega * (csw[ii][jj] + cnw[ii][jj])/2.0;
                                m_css[ii][jj] += omega * (csw[ii][jj] + cse[ii][jj])/2.0;
                    }
                    }
             */
            /*            // 3) According to Reusken
                        REAL alfa = 1.0;
                        for (jj = st; jj <= m_Nx2 + 2; jj += tw)
                        {
                    for (ii = 0; ii <= m_Nx1 + 2; ii += tw)
                            {
                                m_ccc[ii][jj] -= alfa * ( cne[ii][jj] + cnw[ii][jj] +
                                                          cse[ii][jj] + csw[ii][jj]);
                                m_cee[ii][jj] += alfa * (cse[ii][jj] + cne[ii][jj]);
                                m_cnn[ii][jj] += alfa * (cnw[ii][jj] + cne[ii][jj]);
                                m_cww[ii][jj] += alfa * (csw[ii][jj] + cnw[ii][jj]);
                                m_css[ii][jj] += alfa * (csw[ii][jj] + cse[ii][jj]);
                    }
                    }
                    for (jj = 0; jj <= m_Nx2 + 2; jj += tw)
                        {
                    for (ii = st; ii <= m_Nx1 + 2; ii += tw)
                        {
                                m_ccc[ii][jj] -= alfa * ( cne[ii][jj] + cnw[ii][jj] +
                                                          cse[ii][jj] + csw[ii][jj]);
                                m_cee[ii][jj] += alfa * (cse[ii][jj] + cne[ii][jj]);
                                m_cnn[ii][jj] += alfa * (cnw[ii][jj] + cne[ii][jj]);
                                m_cww[ii][jj] += alfa * (csw[ii][jj] + cnw[ii][jj]);
                                m_css[ii][jj] += alfa * (csw[ii][jj] + cse[ii][jj]);
                    }
                    }
             */
        }

        // Voor de overblijvende witte punten hebben we eveneens 
        //                      cnw cnn cne 
        // een 9-punts stencil  cww ccc cee 
        //                      csw css cse 
        // met dus cnw, cne, ccc, csw en cse als witte punten. 
        // Op de witte punten denken we een nieuw schaakbord aangebracht => 
        // S22 en S33. We elimineren nu alle zwarte punten met de zojuist 
        // afgeleide 5-punts formule voor de zwarte punten. We krijgen dan 
        //                          cnn 
        //                      cnw     cne 
        //                  cww     ccc     cee 
        //                      csw     cse 
        //                          css 
        for (jj = st; jj <= m_Nx2 + 2; jj += tw)
        {
            for (ii = st; ii <= m_Nx1 + 2; ii += tw)
            {
                int iw = ii - st;
                int ie = ii + st; //EW: moet dit eerder gedeclareerd worden?
                int js = jj - st;
                int jn = jj + st;
                if (jn <= m_Nx2 + 2) //EW: bijna altijd zo
                {
                    cne[ii][jj] -= m_cnn[ii][jj] * m_cee[ii][jn] / m_ccc[ii][jn];
                    cnw[ii][jj] -= m_cnn[ii][jj] * m_cww[ii][jn] / m_ccc[ii][jn];
                    m_ccc[ii][jj] -= m_cnn[ii][jj] * m_css[ii][jn] / m_ccc[ii][jn];
                    m_cnn[ii][jj] *= -m_cnn[ii][jn] / m_ccc[ii][jn];
                }
                if (js > 0) //EW: alleen bij eerste jj niet zo
                {
                    cse[ii][jj] -= m_css[ii][jj] * m_cee[ii][js] / m_ccc[ii][js];
                    csw[ii][jj] -= m_css[ii][jj] * m_cww[ii][js] / m_ccc[ii][js];
                    m_ccc[ii][jj] -= m_css[ii][jj] * m_cnn[ii][js] / m_ccc[ii][js];
                    m_css[ii][jj] *= -m_css[ii][js] / m_ccc[ii][js];
                }
                if (ie <= m_Nx1 + 2) //EW: bijna altijd zo
                {
                    cne[ii][jj] -= m_cee[ii][jj] * m_cnn[ie][jj] / m_ccc[ie][jj];
                    cse[ii][jj] -= m_cee[ii][jj] * m_css[ie][jj] / m_ccc[ie][jj];
                    m_ccc[ii][jj] -= m_cee[ii][jj] * m_cww[ie][jj] / m_ccc[ie][jj];
                    m_cee[ii][jj] *= -m_cee[ie][jj] / m_ccc[ie][jj];
                }
                if (iw > 0) //EW: alleen bij eerste ii niet zo
                {
                    cnw[ii][jj] -= m_cww[ii][jj] * m_cnn[iw][jj] / m_ccc[iw][jj];
                    csw[ii][jj] -= m_cww[ii][jj] * m_css[iw][jj] / m_ccc[iw][jj];
                    m_ccc[ii][jj] -= m_cww[ii][jj] * m_cee[iw][jj] / m_ccc[iw][jj];
                    m_cww[ii][jj] *= -m_cww[iw][jj] / m_ccc[iw][jj];
                }
            }
        }
        for (jj = 0; jj <= m_Nx2 + 2; jj += tw)
        {
            for (ii = 0; ii <= m_Nx1 + 2; ii += tw)
            {
                int iw = ii - st;
                int ie = ii + st;
                int js = jj - st;
                int jn = jj + st;
                if (jn <= m_Nx2 + 2)
                {
                    cne[ii][jj] -= m_cnn[ii][jj] * m_cee[ii][jn] / m_ccc[ii][jn];
                    cnw[ii][jj] -= m_cnn[ii][jj] * m_cww[ii][jn] / m_ccc[ii][jn];
                    m_ccc[ii][jj] -= m_cnn[ii][jj] * m_css[ii][jn] / m_ccc[ii][jn];
                    m_cnn[ii][jj] *= -m_cnn[ii][jn] / m_ccc[ii][jn];
                }
                if (js > 0)
                {
                    cse[ii][jj] -= m_css[ii][jj] * m_cee[ii][js] / m_ccc[ii][js];
                    csw[ii][jj] -= m_css[ii][jj] * m_cww[ii][js] / m_ccc[ii][js];
                    m_ccc[ii][jj] -= m_css[ii][jj] * m_cnn[ii][js] / m_ccc[ii][js];
                    m_css[ii][jj] *= -m_css[ii][js] / m_ccc[ii][js];
                }
                if (ie <= m_Nx1 + 2)
                {
                    cne[ii][jj] -= m_cee[ii][jj] * m_cnn[ie][jj] / m_ccc[ie][jj];
                    cse[ii][jj] -= m_cee[ii][jj] * m_css[ie][jj] / m_ccc[ie][jj];
                    m_ccc[ii][jj] -= m_cee[ii][jj] * m_cww[ie][jj] / m_ccc[ie][jj];
                    m_cee[ii][jj] *= -m_cee[ie][jj] / m_ccc[ie][jj];
                }
                if (iw > 0)
                {
                    cnw[ii][jj] -= m_cww[ii][jj] * m_cnn[iw][jj] / m_ccc[iw][jj];
                    csw[ii][jj] -= m_cww[ii][jj] * m_css[iw][jj] / m_ccc[iw][jj];
                    m_ccc[ii][jj] -= m_cww[ii][jj] * m_cee[iw][jj] / m_ccc[iw][jj];
                    m_cww[ii][jj] *= -m_cww[iw][jj] / m_ccc[iw][jj];
                }
            }
        }

        // Voor de overblijvende witte punten hebben we 
        //                          cnn 
        //                      cnw     cne 
        //                   cww    ccc    cee 
        //                      csw     cse 
        //                          css 
        // Hierin lumpen we de locale coefficienten cnloc etc. op de rest 
        //  => scheve 5-punts formule 
        //                      cnw     cne 
        //                          ccc 
        //                      csw     cse 

        // Several lumping strategies
        for (jj = st; jj <= m_Nx2 + 2; jj += tw)
        {
            for (ii = st; ii <= m_Nx1 + 2; ii += tw)
            {

                // 1) Similar to RIC
                m_ccc[ii][jj] = m_ccc[ii][jj] + omega * (
                                                         m_cnn[ii][jj] +
                                                         m_cee[ii][jj] +
                                                         m_css[ii][jj] +
                                                         m_cww[ii][jj]);
                /*
                                // 2) Similar to RIC, but on outer elements
                                cne[ii][jj] += omega * (m_cee[ii][jj] + m_cnn[ii][jj])/2.0;
                                cnw[ii][jj] += omega * (m_cww[ii][jj] + m_cnn[ii][jj])/2.0;
                                csw[ii][jj] += omega * (m_cww[ii][jj] + m_css[ii][jj])/2.0;
                                cse[ii][jj] += omega * (m_cee[ii][jj] + m_css[ii][jj])/2.0;
                 */
                /*                // 3) According to Reusken
                                REAL alfa = 1.0;
                        m_ccc[ii][jj] = m_ccc[ii][jj] - alfa * (
                                                m_cnn[ii][jj] +
                                                m_cee[ii][jj] +
                                                m_css[ii][jj] +
                                            m_cww[ii][jj]);
                                cne[ii][jj] += alfa * (m_cee[ii][jj] + m_cnn[ii][jj]);
                                cnw[ii][jj] += alfa * (m_cww[ii][jj] + m_cnn[ii][jj]);
                                csw[ii][jj] += alfa * (m_cww[ii][jj] + m_css[ii][jj]);
                                cse[ii][jj] += alfa * (m_cee[ii][jj] + m_css[ii][jj]);
                 */
                // Om overschrijving te voorkomen gaan we cnw(ii,jj) opslaan in cnn(ii,jj);
                //                                        cne(ii,jj) opslaan in cee(ii,jj);
                //                                        csw(ii,jj) opslaan in cww(ii,jj);
                //                                        cse(ii,jj) opslaan in css(ii,jj);
                m_cnn[ii][jj] = cnw[ii][jj];
                m_cee[ii][jj] = cne[ii][jj];
                m_cww[ii][jj] = csw[ii][jj];
                m_css[ii][jj] = cse[ii][jj];
            }
        }
        ++niv;

        // Uitgaande van het stencil 
        //                          cnn 
        //                     cnw      cne 
        //                 cww      ccc      cee 
        //                     csw      cse 
        //                          css 
        // elimineren we cnw en cne met de scheve 5-punts formule. Dit geeft 
        //                 cnw      cnn      cne 
        //                 cww      ccc      cee 
        //                 csw      css      cse 
        // We zijn nu weer terug de eerste 9-punts formule, maar op een 
        // grover rooster. 
        for (jj = 0; jj <= m_Nx2 + 2; jj += tw)
        {
            for (ii = 0; ii <= m_Nx1 + 2; ii += tw)
            {
                int iw = ii - st;
                int ie = ii + st;
                int js = jj - st;
                int jn = jj + st;
                if (jn <= m_Nx2 + 2 && ie <= m_Nx1 + 2) //EW: bijna altijd zo
                {
                    m_cnn[ii][jj] = m_cnn[ii][jj] - cne[ii][jj] * cnw[ie][jn] / m_ccc[ie][jn];
                    m_cee[ii][jj] = m_cee[ii][jj] - cne[ii][jj] * cse[ie][jn] / m_ccc[ie][jn];
                    m_ccc[ii][jj] = m_ccc[ii][jj] - cne[ii][jj] * csw[ie][jn] / m_ccc[ie][jn];
                    cne[ii][jj] = -cne[ii][jj] * cne[ie][jn] / m_ccc[ie][jn];
                }
                if (iw > 0 && jn <= m_Nx2 + 2) //EW: bijna altijd zo
                {
                    m_cnn[ii][jj] = m_cnn[ii][jj] - cnw[ii][jj] * cne[iw][jn] / m_ccc[iw][jn];
                    m_cww[ii][jj] = m_cww[ii][jj] - cnw[ii][jj] * csw[iw][jn] / m_ccc[iw][jn];
                    m_ccc[ii][jj] = m_ccc[ii][jj] - cnw[ii][jj] * cse[iw][jn] / m_ccc[iw][jn];
                    cnw[ii][jj] = -cnw[ii][jj] * cnw[iw][jn] / m_ccc[iw][jn];
                }
                if (js > 0 && iw > 0) //EW: bijna altijd zo
                {
                    m_css[ii][jj] = m_css[ii][jj] - csw[ii][jj] * cse[iw][js] / m_ccc[iw][js];
                    m_cww[ii][jj] = m_cww[ii][jj] - csw[ii][jj] * cnw[iw][js] / m_ccc[iw][js];
                    m_ccc[ii][jj] = m_ccc[ii][jj] - csw[ii][jj] * cne[iw][js] / m_ccc[iw][js];
                    csw[ii][jj] = -csw[ii][jj] * csw[iw][js] / m_ccc[iw][js];
                }
                if (js > 0 && ie <= m_Nx1 + 2) //EW: bijna altijd zo
                {
                    m_css[ii][jj] = m_css[ii][jj] - cse[ii][jj] * csw[ie][js] / m_ccc[ie][js];
                    m_cee[ii][jj] = m_cee[ii][jj] - cse[ii][jj] * cne[ie][js] / m_ccc[ie][js];
                    m_ccc[ii][jj] = m_ccc[ii][jj] - cse[ii][jj] * cnw[ie][js] / m_ccc[ie][js];
                    cse[ii][jj] = -cse[ii][jj] * cse[ie][js] / m_ccc[ie][js];
                }
            }
        }
        st = tw;
        tw *= 2;

    } // end of RRB-iteration

    /* The last Schur complement will be solved exactly with a complete Cholesky decomposition.
     * This decomposition Smax=LL^T will be stored in Chol=L.
     * The matrix Chol will have dimension of the bandwith*schurdimension,
     * in case of a five point stencil, the bandwith is NSx1 or NSx2, with NSx1 and NSx2 the dimensions of
     * the coarsest grid, so schurdimension = NSx1 * NSx2.
     * The coarsest grid is given by the last loop of the RRB-iteration. Be carefull, because after
     * the loop st=tw and tw*=2 is done.
     * So NSx1 = floor( (Nx1+3)/st ) and NSx2 = floor( (Nx2+3)/st ). */

    /* Calculate the dimensions of the coarsest grid:
     * NSx1 = floor( (Nx1+3)/st ) and NSx2 = floor( (Nx2+3)/st ).
     * Dimension of the matrix will be NSx1*NSx2 x NSx1*NSx2 = dimS x dimS.
     * The bandwith depends on horizontal/vertical numbering: either NSx1 or NSx2.
     * The best choice is min(NSx1,NSx2), but take now NSx1 (=horizontal numbering).
     * The diagonals of the complete Cholesky decomposition are stored in U, with
     * dimensions dimS x NSx1+1. */

    // Mesh size dN of coarsest grid.
    // EW: dit is eigenlijk niet nodig, om dN uit te rekenen, gebruik st.
    //    int dN = pow(2.0,m_nivmax-1);

    // Number of nodes of the coarsest grid.
    NSx1 = (m_Nx1 + 2) / st + 1;
    NSx2 = (m_Nx2 + 2) / st + 1;
    m_dimS = NSx1 * NSx2;

    // The bandwith of the matrix on the coarsest grid.
    /* EW: hier moet eigenlijk min(NSx1,NSx2) gelden. Maar dan moet ook de rest meeveranderd worden.
     * Dit komt neer op de keuze horizontale/verticale nummering. */
    m_bandwith = NSx1;

    // For the Cholesky decomposition. Only the diagonals are stored.
    m_Chol.resize(m_dimS, m_bandwith + 1);

    /* Store the diagonals in Chol: Chol(:,0) is the main diagonal, Chol(:,k) is the k-th left diagonal.
     * Also lump de nine-point stencil to a five-point stencil. */
#ifndef EWNINEPOINT
    int kk = 0; // counting variable
    for (jj = 0; jj <= m_Nx2 + 2; jj += st)
    {
        for (ii = 0; ii <= m_Nx1 + 2; ii += st)
        {
            m_Chol[kk][0] = m_ccc[ii][jj] + omega * (cne[ii][jj] + cnw[ii][jj]
                                                     + csw[ii][jj] + cse[ii][jj]);
            // EW: dit vorige hoort eigenlijk tot de RRB-procedure

            m_Chol[kk][1] = m_cww[ii][jj];
            m_Chol[kk][m_bandwith] = m_css[ii][jj];

            // For the easy implementation of an upcoming loop, set ccc to one.
            m_ccc[ii][jj] = 1.0;

            ++kk;
        }
    }
#else
    m_Chol.resize(m_dimS, m_bandwith + 2);
    int kk = 0; // counting variable
    for (jj = 0; jj <= m_Nx2 + 2; jj += st)
    {
        for (ii = 0; ii <= m_Nx1 + 2; ii += st)
        {
            m_Chol[kk][0] = m_ccc[ii][jj];
            // EW: dit vorige hoort eigenlijk tot de RRB-procedure

            m_Chol[kk][1] = m_cww[ii][jj];
            m_Chol[kk][m_bandwith - 1] = cse[ii][jj];
            m_Chol[kk][m_bandwith] = m_css[ii][jj];
            m_Chol[kk][m_bandwith + 1] = csw[ii][jj];

            // For the easy implementation of an upcoming loop, set ccc to one.
            m_ccc[ii][jj] = 1.0;

            ++kk;
        }
    }
#endif

#ifndef EWLAPACK
    // Perform the complete Cholesky decomposition of the last part.
    // EW: zie Cholesky.m; let op dat in Matlab matrices bij (1,1) beginnen, hier bij (0,0).
    sumc = 0.0; // dummy variable for summation during Cholesky decomposition
    for (jj = 0; jj < m_dimS; ++jj)
    {
        sumc = 0.0;
        for (kk = std::max(1, jj + 1 - m_bandwith); kk <= jj; ++kk)
        {
            sumc += FloatUtils::square(m_Chol[jj][jj + 1 - kk]);
        }
        m_Chol[jj][0] = sqrt(m_Chol[jj][0] - sumc);
        int iiend = std::min(jj + 1 + m_bandwith, m_dimS); // EW: is dit echt sneller?
        for (ii = jj + 1; ii < iiend; ++ii)
        {
            sumc = 0.0;
            int kkbegin = std::max(1, ii + 1 - m_bandwith); // EW: wellicht eerder declareren
            for (kk = kkbegin; kk <= jj; ++kk)
            {
                sumc += m_Chol[ii][ii + 1 - kk] * m_Chol[jj][jj + 1 - kk];
            }
            m_Chol[ii][ii - jj] = (m_Chol[ii][ii - jj] - sumc) / m_Chol[jj][0];
        }
    }

#else

    // Resize and initialize the matrix m_LaChol.
    // EW: Wellicht bwE*2+1 als bandbreedte.
#ifndef EWNINEPOINT
    m_LaChol.resize(m_dimS, m_bandwith);
#else
    m_LaChol.resize(m_dimS, m_bandwith + 1);
#endif
    m_LaChol = 0.0;

    // Subsitute the matrix m_Chol (Array2D<REAL>) to matrix m_LaChol (LaSymmBandMatdouble).
    for (ii = 0; ii < m_dimS; ++ii)
    {
        m_LaChol(ii, ii) = m_Chol[ii][0]; // center elements on main diagonal
    }
    for (ii = 1; ii < m_dimS; ++ii)
    {
        m_LaChol(ii, ii - 1) = m_Chol[ii][1]; // west elements on first outer diagonal
    }
    for (ii = m_bandwith; ii < m_dimS; ++ii)
    {
        m_LaChol(ii, ii - m_bandwith) = m_Chol[ii][m_bandwith]; // south elements on bwE-th outer diagonal
    }
#ifdef EWNINEPOINT
    for (ii = m_bandwith - 1; ii < m_dimS; ++ii)
    {
        m_LaChol(ii, ii - m_bandwith + 1) = m_Chol[ii][m_bandwith - 1]; // southeast elements on bwE-1-th outer diagonal
    }
    for (ii = m_bandwith + 1; ii < m_dimS; ++ii)
    {
        m_LaChol(ii, ii - m_bandwith - 1) = m_Chol[ii][m_bandwith + 1]; // southwest elements on bwE+1-th outer diagonal
    }
#endif

    // Calculate the Cholesky decomposition.
    LaSymmBandMatFactorizeIP(m_LaChol);
#endif


    // Achteraf delen door hoofddiagonaalelement om 
    // incomplete LDL'-ontbinding te krijgen:
    // EW: moet dit echt over alle punten lopen?
    for (ii = 0; ii <= m_Nx1 + 2; ++ii)
    {
        for (jj = 0; jj <= m_Nx2 + 2; ++jj)
        {
            m_ccc[ii][jj] = 1.0 / m_ccc[ii][jj];
            m_cnn[ii][jj] = m_cnn[ii][jj] * m_ccc[ii][jj];
            m_css[ii][jj] = m_css[ii][jj] * m_ccc[ii][jj];
            m_cee[ii][jj] = m_cee[ii][jj] * m_ccc[ii][jj];
            m_cww[ii][jj] = m_cww[ii][jj] * m_ccc[ii][jj];
        }
    }

    // Om deling door nul te voorkomen wordt voor alle witte punten 
    // van een schaakbordvolgorde niet de hoofddiagonaal van 
    // de matrix opgeslagen, maar 1/(--). 
    /* The center stencil will be used in the matvec() as inverse for the black points and
     * as itself in the red points of first level. Therefore, the center stencil is stored
     * as its inverse only on the black points.
     * Also in the elimination of all first order black points, the inverse is needed. */
    for (jj = 0; jj <= m_Nx2 + 2; jj += 2)
    {
        for (ii = 1; ii <= m_Nx1 + 2; ii += 2)
        {
            cc[ii][jj] = 1.0 / cc[ii][jj];
        }
    }
    for (jj = 1; jj <= m_Nx2 + 2; jj += 2)
    {
        for (ii = 0; ii <= m_Nx1 + 2; ii += 2)
        {
            cc[ii][jj] = 1.0 / cc[ii][jj];
        }
    }
    // De complete Choleski-ontbinding van het laatste stuk is nu gemaakt.

} // end of preconditioner


/// solves preconditioner

/* During the CG-iteration the equation z = M^-1 r has to be calculated, with M the
 * preconditioner. Because of the RRB-structure of the preconditioner, an LDL^T decomposition
 * is performed. So z = (LDL^T)^-1 r is solved with forward and backward substitution.
 * The complete Cholesky decomposition at the coarsest grid will be solved and the solution
 * of it is substituted in the D of the decomposition. */
void RrbElwinSolver::precsolve(Array2D<REAL> *pz)
{
    Array2D<REAL> &z = *pz; // input-vector as well as output-vector: z = M^-1 z
    Array1D<REAL> rhs(m_dimS); // right hand side of complete Cholesky decomposition
    Array1D<REAL> sol(m_dimS); // solution of complete Cholesky decomposition
    // EW: deze twee variabelen zijn niet strikt noodzakelijk
    // EW: TO DO: dit kunnen ook Array1D<REAL> zijn.

    int ii, jj; // counting variables

    int st, tw, niv; // dummy variables during the forward and backward substitution
    // Eerst oplossen L(Az)=Az. We maken hierbij gebruik van symmetrie.
    // EW: voor exact staartdeel: niv < nivmax
    st = 1;
    tw = 2;
    niv = 1;
    while (niv < m_nivmax)
    {
        if (niv > 1)
        {
            // EW: als het laatste rooster maar 'e'en punt heeft, krijg je hier problemen.
            // Do the forward substitution for the points (0,jj): black points of level niv at west boundary
            ii = 0;
            for (jj = st; jj <= m_Nx2 + 1; jj += tw)
            {
                z[ii ][jj - st] -= m_css[ii][jj] * z[ii][jj];
                z[ii + st][jj ] -= m_cee[ii][jj] * z[ii][jj];
                z[ii ][jj + st] -= m_cnn[ii][jj] * z[ii][jj];
            }

            // Do the forward substitution for the points (ii,0): black points of level niv at south boundary
            jj = 0;
            for (ii = st; ii <= m_Nx1 + 1; ii += tw)
            {
                z[ii + st][jj ] -= m_cee[ii][jj] * z[ii][jj];
                z[ii ][jj + st] -= m_cnn[ii][jj] * z[ii][jj];
                z[ii - st][jj ] -= m_cww[ii][jj] * z[ii][jj];
            }

            // Do the forward substitution for the points (ii,jj): black points of level niv outside boundaries
            for (jj = st; jj <= m_Nx2 + 1; jj += tw)
            {
                for (ii = tw; ii <= m_Nx1 + 1; ii += tw)
                {
                    z[ii ][jj - st] -= m_css[ii][jj] * z[ii][jj];
                    z[ii + st][jj ] -= m_cee[ii][jj] * z[ii][jj];
                    z[ii ][jj + st] -= m_cnn[ii][jj] * z[ii][jj];
                    z[ii - st][jj ] -= m_cww[ii][jj] * z[ii][jj];
                }
            }
            for (jj = tw; jj <= m_Nx2 + 1; jj += tw)
            {
                for (ii = st; ii <= m_Nx1 + 1; ii += tw)
                {
                    z[ii ][jj - st] -= m_css[ii][jj] * z[ii][jj];
                    z[ii + st][jj ] -= m_cee[ii][jj] * z[ii][jj];
                    z[ii ][jj + st] -= m_cnn[ii][jj] * z[ii][jj];
                    z[ii - st][jj ] -= m_cww[ii][jj] * z[ii][jj];
                }
            }
        }

        // Do the forward substitution for the points (ii,jj): red points of level niv not in level niv+1
        for (jj = st; jj <= m_Nx2 + 1 - st; jj += tw)
        {
            for (ii = st; ii <= m_Nx1 + 1 - st; ii += tw)
            {
                z[ii + st][jj - st] -= m_css[ii][jj] * z[ii][jj];
                z[ii - st][jj - st] -= m_cww[ii][jj] * z[ii][jj];
                z[ii + st][jj + st] -= m_cee[ii][jj] * z[ii][jj];
                z[ii - st][jj + st] -= m_cnn[ii][jj] * z[ii][jj];
            }
        }

        ++niv; //niv += 2;
        st = tw;
        tw *= 2;
    } // end of L(Az)=Az

    // Nu oplossen D(Az)=Az:
    for (jj = 1; jj <= m_Nx2 + 2; jj += 2)
    {
        for (ii = 1; ii <= m_Nx1 + 2; ii += 2)
        {
            z[ii ][jj ] *= m_ccc[ii ][jj ];
            z[ii - 1][jj - 1] *= m_ccc[ii - 1][jj - 1];
        }
    }

    // At the coarsest grid, solve the Cholesky decomposition.
    // The right hand side
    rhs = 0.0;

    /* Fill the right hand side of the LDL^T decomposition in rhs, the
     * right hand side of the Cholesky decomposition */

    // Mesh size dN of coarsest grid.
    // EW: dit is eigenlijk niet nodig, om dN uit te rekenen, gebruik st.
    //    int dN = pow(2.0,m_nivmax-1);

#ifndef EWLAPACK
    int kk = 0;
    for (jj = 0; jj <= m_Nx2 + 2; jj += st)
    {
        for (ii = 0; ii <= m_Nx1 + 2; ii += st)
        {
            rhs[kk] = z[ii][jj];
            ++kk;
        }
    }

#else

    LaVectorDouble v(m_dimS);
    int kk = 0;
    for (jj = 0; jj <= m_Nx2 + 2; jj += st)
    {
        for (ii = 0; ii <= m_Nx1 + 2; ii += st)
        {
            v(kk) = z[ii][jj];
            ++kk;
        }
    }
#endif

#ifndef EWLAPACK
    // The solution vector
    sol = 0.0;
    // EW: in principe kan sol over rhs geschreven worden.
    // forward substitution to solve L y = b.
    for (ii = 0; ii < m_dimS; ++ii)
    {
        sol[ii] = rhs[ii];
        for (jj = std::max(0, ii - m_bandwith); jj < ii; ++jj)
        {
            sol[ii] -= m_Chol[ii][ii - jj] * sol[jj];
        }
        sol[ii] /= m_Chol[ii][0];
    }
    // backward substitution to solve L^T x = y.
    for (ii = m_dimS - 1; ii >= 0; --ii)
    {
        for (jj = ii + 1; jj <= std::min(m_dimS - 1, ii + m_bandwith); ++jj)
        {
            sol[ii] -= m_Chol[jj][jj - ii] * sol[jj];
        }
        sol[ii] /= m_Chol[ii][0];
    }

#else

    // Solve the system.
    LaLinearSolveIPChol(m_LaChol, v);
#endif

    /* Substitute the solution of the Cholesky decomposition in 
     * the solution vector of the LDL^T decomposition. */
#ifndef EWLAPACK
    kk = 0;
    for (jj = 0; jj <= m_Nx2 + 2; jj += st)
    {
        for (ii = 0; ii <= m_Nx1 + 2; ii += st)
        {
            z[ii][jj] = sol[kk];
            ++kk;
        }
    }

#else

    kk = 0;
    for (jj = 0; jj <= m_Nx2 + 2; jj += st)
    {
        for (ii = 0; ii <= m_Nx1 + 2; ii += st)
        {
            z[ii][jj] = v(kk);
            ++kk;
        }
    }
#endif


    // Nu oplossen U(Az)=Az:
    while (niv > 1)
    {
        tw = st; // The first 'st' originates from the L(Az)=Az loop.
        st /= 2;
        --niv; //niv += -2; //EW niv -= 2

        // Do the backward substitution for the points (ii,jj): red points of level niv not in level niv+1
        for (jj = st; jj <= m_Nx2 + 1; jj += tw)
        {
            for (ii = st; ii <= m_Nx1 + 1; ii += tw)
            {
                z[ii][jj] += -m_css[ii][jj] * z[ii + st][jj - st]
                    - m_cww[ii][jj] * z[ii - st][jj - st]
                    - m_cee[ii][jj] * z[ii + st][jj + st]
                    - m_cnn[ii][jj] * z[ii - st][jj + st];
            }
        }

        // Do the backward substitution for the black points of level niv
        if (niv > 1)
        {
            // EW: als het laatste rooster maar 'e'en punt heeft, krijg je hier problemen.
            // Do the backward substitution for the points (0,jj): black points of level niv at west boundary
            ii = 0;
            for (jj = st; jj <= m_Nx2 + 1; jj += tw)
            {
                z[ii][jj] += -m_css[ii][jj] * z[ii ][jj - st]
                    - m_cee[ii][jj] * z[ii + st][jj ]
                    - m_cnn[ii][jj] * z[ii ][jj + st];
            }

            // Do the backward substitution for the points (ii,0): black points of level niv at south boundary
            jj = 0;
            for (ii = st; ii <= m_Nx1 + 1; ii += tw)
            {
                z[ii][jj] += -m_cww[ii][jj] * z[ii - st][jj ]
                    - m_cee[ii][jj] * z[ii + st][jj ]
                    - m_cnn[ii][jj] * z[ii ][jj + st];
            }

            // Do the backward substitution for the points (ii,jj): black points of level niv outside boundaries
            for (jj = st; jj <= m_Nx2 + 1; jj += tw)
            {
                for (ii = tw; ii <= m_Nx1 + 1; ii += tw)
                {
                    z[ii][jj] += -m_css[ii][jj] * z[ii ][jj - st]
                        - m_cww[ii][jj] * z[ii - st][jj ]
                        - m_cee[ii][jj] * z[ii + st][jj ]
                        - m_cnn[ii][jj] * z[ii ][jj + st];
                }
            }
            for (jj = tw; jj <= m_Nx2 + 1; jj += tw)
            {
                for (ii = st; ii <= m_Nx1 + 1; ii += tw)
                {
                    z[ii][jj] += -m_css[ii][jj] * z[ii ][jj - st]
                        - m_cww[ii][jj] * z[ii - st][jj ]
                        - m_cee[ii][jj] * z[ii + st][jj ]
                        - m_cnn[ii][jj] * z[ii ][jj + st];
                }
            }
        }
    } // end of U(Az)=Az


} // end of precsolve


/// solves system
/*  In deze subroutine wordt S(x)=R opgelost, waarbij S het 
    Schurcomplement is dat ontstaat door alle 
    zwarte punten van een schaakbordvolgorde te elimineren. 
    We stoppen de iteratie zodra de 2-norm van het 
    residu tenminste een factor tol is afgenomen. 
    Het CG-algoritme is geimplementeerd zoals op blz. 
    146 in 'Solving Linear Systems on Vector and Shared 
    Memory Computers' van Dongarra, Duff, Sorensen en 
    v.d. Vorst is beschreven. 
    parameters:
    X: veld waarin de oplossing komt te staan. 
    B: veld met rechterlid van het stelsel vergelijkingen. 
       Deze vector wordt overschreven door het residu.  */
// EW: wordt niet het hele systeem opgelost?

void RrbElwinSolver::solve(Array2D<REAL> *pB, // real B(0:Nx1+2,0:Nx2+2)
                           Array2D<REAL> *pX) // real X(0:Nx1+2,0:Nx2+2)
{
#ifdef _OPENMP
#pragma omp master
#endif
    { //------------------------------------------------------------------------

    Array2D<REAL> &r = *pB; // right hand side linear system and residual during CG-iteration
    Array2D<REAL> &x = *pX; // initial solution linear system and estimate during CG-iteration

    Array2D<REAL> &cc = *m_pcc; // center stencil linear system
    Array2D<REAL> &cs = *m_pcs; // south stencil linear system
    Array2D<REAL> &cw = *m_pcw; // west stencil linear system

    Array2D<REAL> p(r.dim1(), r.dim2()); // search direction
    Array2D<REAL> q(r.dim1(), r.dim2()); // multiplied search vector
    Array2D<REAL> y(r.dim1(), r.dim2()); // dummy for matrix-vector multiplication

    /* In the loops of precsolve(), the indices will vary according to the RRB-orderings. To prevent
     * if-statements at the north and east boundaries, the array z will be enlarged.
     * It has to have dimension r.dim1+st_max and r.dim1+st_max with st_max = 2^(nivmax-1). */
    /* EW: the solver could be better (less memory) when z is not enlarged, at the expense of
     * some nasty if-statements in precsolve(). */
    int stmax = pow(2.0, m_nivmax - 1);
    Array2D<REAL> z(r.dim1() + stmax, r.dim2() + stmax); // preconditioned residual

    REAL stopvalue, alpha, beta, rhoold, rhonew, sigma;

    int ii, jj, iteration; // counting variables

    /* EW: In deze routine wordt ook cne, cnw en cns gebruikt. Dit zijn niet de stencils,
     * maar hier staan variabelen voor CG in, zoals p en q, de zoekrichting. Dus erg
     * verwarrend en waarschijnlijk onnodig (gebruik lokale variabelen). */

    /* EW: Voor i,j=1,3,5,... wordt het eerste Schur-complement opgelost. De andere helft
     * punten worden met een Gauss-eliminatie slag exact uitgerekend. */

    // Alle zwarte punten worden exact geelimineerd:
    /* Calculate the right hand side of the first Schur complement. This is done by eliminating
     * all black points of first level: b=b-Sb. */

    // Do this elimination for all first level red points in the domain without ghost points.
    for (jj = 1; jj <= m_Nx2 + 1; jj += 2)
    {
        for (ii = 1; ii <= m_Nx1 + 1; ii += 2)
        {
            r[ii][jj] += -cs[ii ][jj + 1] * r[ii ][jj + 1] * cc[ii ][jj + 1]
                - cw[ii + 1][jj ] * r[ii + 1][jj ] * cc[ii + 1][jj ]
                - cs[ii ][jj ] * r[ii ][jj - 1] * cc[ii ][jj - 1]
                - cw[ii ][jj ] * r[ii - 1][jj ] * cc[ii - 1][jj ];
        }
    }
    for (jj = 2; jj <= m_Nx2 + 1; jj += 2)
    {
        for (ii = 2; ii <= m_Nx1 + 1; ii += 2)
        {
            r[ii][jj] += -cs[ii ][jj + 1] * r[ii ][jj + 1] * cc[ii ][jj + 1]
                - cw[ii + 1][jj ] * r[ii + 1][jj ] * cc[ii + 1][jj ]
                - cs[ii ][jj ] * r[ii ][jj - 1] * cc[ii ][jj - 1]
                - cw[ii ][jj ] * r[ii - 1][jj ] * cc[ii - 1][jj ];
        }
    }
    // Do this elimination for point (0,0): southwest corner of first level red points.
    ii = 0;
    jj = 0;
    r[ii][jj] += -cs[ii ][jj + 1] * r[ii ][jj + 1] * cc[ii ][jj + 1]
        - cw[ii + 1][jj ] * r[ii + 1][jj ] * cc[ii + 1][jj ];
    // Do this elimination for points (ii,0): south boundary of first level red points.
    for (ii = 2; ii <= m_Nx1 + 1; ii += 2)
    {
        r[ii][jj] += -cs[ii ][jj + 1] * r[ii ][jj + 1] * cc[ii ][jj + 1]
            - cw[ii + 1][jj ] * r[ii + 1][jj ] * cc[ii + 1][jj ]
            - cw[ii ][jj ] * r[ii - 1][jj ] * cc[ii - 1][jj ];
    }
    // Do this elimination for points (0,jj): west boundary of first level red points.
    ii = 0;
    for (jj = 2; jj <= m_Nx2 + 1; jj += 2)
    {
        r[ii][jj] += -cs[ii ][jj + 1] * r[ii ][jj + 1] * cc[ii ][jj + 1]
            - cw[ii + 1][jj ] * r[ii + 1][jj ] * cc[ii + 1][jj ]
            - cs[ii ][jj ] * r[ii ][jj - 1] * cc[ii ][jj - 1];
    }
    // The black points are eliminated

    // We berekenen de vector S*X met S het 1e Schurcomplement.
    /* Input: stencils cc, cs and cw, dimension Nx1 and Nx2, and initial solution x.
     * Output: arrays m_cse=S*x=Q2*psi0 and m_cne */
    // EW: hier wordt over alle punten iets uitgerekend. CG over de helft ervan.
    //matvec(cc, cs, cw, x, &m_cse, &m_cne, m_Nx1, m_Nx2);

    // Calculate q=Sx. Now q is not the preconditioned search vector, but S*x.
    //y = 0.0;
    matvec(x, &q, &y);

    // Calculate residual: r=b-A*x: b=b-S*x.
    /* EW: Kan je cne=0.0 gebruiken? Er wordt hier over de helft van de punten geloopt,
     * de rest van de punten in cne komen uit vorige matvec().
     * cne=0 is nodig voor eerste iteratie CG (p^0=0) */
    /* EW: dit soort loops zijn waarschijnlijk de reden voor de tweede rand spookpunten.
     * Dus als dit wordt verbeterd, dan kan deze rand weg. */
    for (jj = 1; jj <= m_Nx2 + 2; jj += 2)
    {
        for (ii = 1; ii <= m_Nx1 + 2; ii += 2)
        {
            // Calculate residual r=b-Ax on all first level red points.
            r[ii - 1][jj - 1] -= q[ii - 1][jj - 1];
            r[ii ][jj ] -= q[ii ][jj ];
        }
    }

    // Solve z = M^-1 r, the preconditioned residual
    for (jj = 0; jj <= m_Nx2 + 2; ++jj)
    {
        for (ii = 0; ii <= m_Nx1 + 2; ++ii)
        {
            z[ii][jj] = r[ii][jj];
        }
    }
    //z = r;
    precsolve(&z);

    // Initialize parameters and calculate the norm of the preconditioned residual: rhonew=<r,z>.
    beta = 0.0;
    rhoold = 0.0;
    rhonew = 0.0;
    for (jj = 1; jj <= m_Nx2 + 2; jj += 2)
    {
        for (ii = 1; ii <= m_Nx1 + 2; ii += 2)
        {
            rhonew += r[ii - 1][jj - 1] * z[ii - 1][jj - 1]
                + r[ii ][jj ] * z[ii ][jj ];
        }
    }

    // The termination criterium for the CG-iteration.
    //stopvalue = FloatUtils::square(m_tolerance);
    //stopvalue = (rhonew + 1) * FloatUtils::square(m_tolerance);
    stopvalue = FloatUtils::square(m_tolerance);
    assert(stopvalue > 0);

    //    stopvalue = FloatUtils::square(2.0e-6);
    // gekozen waarde levert goede resultaten maar kan wellicht nog beter

    // Start the CG-iteration.
    // EW: fabs moet weg kunnen
    for (iteration = 1; fabs(rhonew) > stopvalue && iteration < maxIterations; ++iteration)
    {

        // EW: in eerste iteratie geldt al beta=0. Toch kijken wat er dan in cse en cne zit.
        // Calculating the search vector: p = z + beta p. m_cne contains the search direction.
        // m_cse should contain the residual.
        /* We have p = z + beta * cne, with cne from the matvec() and later on set to zero during
         * the initialization. Also beta=0 in initialization. */ //EW: niet erg netjes gedaan.
        // During the next iterations, cne is the search vector.

        // Calculate the search vector m_cne:
        if (iteration == 1)
        {
            for (jj = 1; jj <= m_Nx2 + 2; jj += 2)
            {
                for (ii = 1; ii <= m_Nx1 + 2; ii += 2)
                {
                    p[ii - 1][jj - 1] = z[ii - 1][jj - 1];
                    p[ii ][jj ] = z[ii ][jj ];
                }
            }
        }
        else
        {
            beta = rhonew / rhoold;
            for (jj = 1; jj <= m_Nx2 + 2; jj += 2)
            {
                for (ii = 1; ii <= m_Nx1 + 2; ii += 2)
                {
                    p[ii - 1][jj - 1] = z[ii - 1][jj - 1] + beta * p[ii - 1][jj - 1];
                    p[ii ][jj ] = z[ii ][jj ] + beta * p[ii ][jj ];
                }
            }
        }

        // We berekenen de vector S*P met S het 1e Schurcomplement en P de zoekrichting: 
        // cc,cs,cw the stencil of the matrix, cne input and cse=Az=q, cnw=y.
        // Calculate q=Sp, with p search vector and S first Schur complement.
        //matvec(cc, cs, cw, m_cne, &m_cse, &m_cnw, m_Nx1, m_Nx2);
        matvec(p, &q, &y);

        // sigma = <p,q>
        sigma = 0.0;
        for (jj = 1; jj <= m_Nx2 + 2; jj += 2)
        {
            for (ii = 1; ii <= m_Nx1 + 2; ii += 2)
            {
                sigma += p[ii - 1][jj - 1] * q[ii - 1][jj - 1]
                    + p[ii ][jj ] * q[ii ][jj ];
            }
        }

        // alpha = <r,r>/<p,q> = rhonew/sigma
        alpha = rhonew / sigma;

        /* Update solution and residual: x = x + alpha p and r = r - alpha q on 
         * all first level red points. Calculate rho = <r,r> the unpreconditioned residual norm. */
        for (jj = 1; jj <= m_Nx2 + 2; jj += 2)
        {
            for (ii = 1; ii <= m_Nx1 + 2; ii += 2)
            {
                x[ii - 1][jj - 1] += alpha * p[ii - 1][jj - 1];
                x[ii ][jj ] += alpha * p[ii ][jj ];

                r[ii - 1][jj - 1] -= alpha * q[ii - 1][jj - 1];
                r[ii ][jj ] -= alpha * q[ii ][jj ];
            }
        }

        z = 0.0; // EW: weet niet zeker of nodig is.
        for (jj = 0; jj <= m_Nx2 + 2; ++jj)
        {
            for (ii = 0; ii <= m_Nx1 + 2; ++ii)
            {
                z[ii][jj] = r[ii][jj];
            }
        }
        //z = r;
        precsolve(&z);

        rhoold = rhonew;
        // rhonew = <r,r>
        rhonew = 0.0;
        for (jj = 1; jj <= m_Nx2 + 2; jj += 2)
        {
            for (ii = 1; ii <= m_Nx1 + 2; ii += 2)
            {
                rhonew += r[ii - 1][jj - 1] * z[ii - 1][jj - 1]
                    + r[ii ][jj ] * z[ii ][jj ];
            }
        }
    } // end of CG-iteration

    // Alle componenten van de eerste groep berekenen door 
    // terugsubstitutie: x1=D^(-1)*(RL1-A12(x2))
    /* The eliminated black nodes are calculated with x1=D1^(-1)*(b1-A*x2), with x2 the
     * nodes calculated in the CG-iteration and x1 the eliminated nodes and b1 the
     * right hand side of the original system and D1 its diagonal.
     * Note that r[ii][jj] in the next loops still contains the original right
     * hand side, cs,cw,cc the original stencil and x the solution from CG.
     * The array cc contains the inverse of the diagonal of the original system.
     * The system is not diagonally scaled.*/

    // Do this substitution for points (ii,jj): all first level black points outside boundaries.
    for (jj = 1; jj <= m_Nx2 + 1; jj += 2)
    {
        for (ii = 2; ii <= m_Nx1 + 1; ii += 2)
        {
            x[ii][jj] = (r[ii][jj]
                         - cs[ii ][jj + 1] * x[ii ][jj + 1]
                         - cw[ii + 1][jj ] * x[ii + 1][jj ]
                         - cw[ii ][jj ] * x[ii - 1][jj ]
                         - cs[ii ][jj ] * x[ii ][jj - 1]) * cc[ii][jj];
        }
    }
    for (jj = 2; jj <= m_Nx2 + 1; jj += 2)
    {
        for (ii = 1; ii <= m_Nx1 + 1; ii += 2)
        {
            x[ii][jj] = (r[ii][jj]
                         - cs[ii ][jj + 1] * x[ii ][jj + 1]
                         - cw[ii + 1][jj ] * x[ii + 1][jj ]
                         - cw[ii ][jj ] * x[ii - 1][jj ]
                         - cs[ii ][jj ] * x[ii ][jj - 1]) * cc[ii][jj];
        }
    }
    // Do this substitution for points (0,jj): south boundary of first level black points.
    ii = 0;
    for (jj = 1; jj <= m_Nx2 + 1; jj += 2)
    {
        x[ii][jj] = (r[ii][jj]
                     - cs[ii ][jj + 1] * x[ii ][jj + 1]
                     - cw[ii + 1][jj ] * x[ii + 1][jj ]
                     - cs[ii ][jj ] * x[ii ][jj - 1]) * cc[ii][jj];
    }
    // Do this substitution for points (ii,0): west boundary of first level black points.
    jj = 0;
    for (ii = 1; ii <= m_Nx1 + 1; ii += 2)
    {
        x[ii][jj] = (r[ii][jj]
                     - cs[ii ][jj + 1] * x[ii ][jj + 1]
                     - cw[ii + 1][jj ] * x[ii + 1][jj ]
                     - cw[ii ][jj ] * x[ii - 1][jj ]) * cc[ii][jj];
    }

    } //master--------------------------------------------------------------
#ifdef _OPENMP
#pragma omp barrier
#ifdef _OPENMP
} // end of solve 


/// matrix vector product
/* Deze subroutine berekent de vector Az=S*x met S het Schurcomplement dat 
 * ontstaat door alle zwarte punten van een schaakbordvolgorde te elimineren.  */
/* Subroutine to compute the matrix-vector multiplication q=Sp, with S the first Schur complement.
 * Because the Schur complement has the form S = D - S D^-1 S, first y = D^-1 S p is calculated.
 * Then q = D p - S y will be calculated.
 * Input:
 *   Member variables: center, south and west stencil (cc,cs,cw) of matrix and size of domain (Nx1,Nx2).
 *   From function call: vector (x) which has to be multiplied by the first Schur complement.
 * Output: vectors q and resulting from the multiplication:
 *   y: the vector y = D^-1 S p on all first order black points;
 *   q: the vector q = D p - S y on all first order red points.
 * Note that the center stencil cc is stored as 1/cc inside the preconditioner(). */
// EW: waarom is het stencil const& en de vectoren een pointer? const variabelen mag je niet veranderen
/* EW: wellicht kan y een lokale variabele worden, in principe wordt y nooit gebruikt. In de CG-iteratie
 * wordt y niet (meer) gebruikt. Toch geen lokale variabele van matvec, want anders moet tijdens elke
 * CG-iteratie de array gedeclareerd worden, nu maar 'e'e keer per solve-procedure. */
// EW: Waarschijnlijk kan q in dezelfde array als y worden gezet.

void RrbElwinSolver::matvec(const Array2D<REAL> &p,
                            Array2D<REAL> *pq,
                            Array2D<REAL> *py)
{

    Array2D<REAL> &q = *pq;
    Array2D<REAL> &y = *py;

    Array2D<REAL> &cc = *m_pcc; // center stencil linear system
    Array2D<REAL> &cs = *m_pcs; // south stencil linear system
    Array2D<REAL> &cw = *m_pcw; // west stencil linear system

    int ii, jj;

    // Calcalate the matrix-vector multiplication y = D^-1 S p.

    // Do this multiplication for points (ii,jj): all first level black points outside boundaries.
    for (jj = 1; jj <= m_Nx2 + 1; jj += 2)
    {
        for (ii = 2; ii <= m_Nx1 + 1; ii += 2)
        {
            y[ii][jj] = (cs[ii ][jj + 1] * p[ii ][jj + 1]
                         + cw[ii + 1][jj ] * p[ii + 1][jj ]
                         + cs[ii ][jj ] * p[ii ][jj - 1]
                         + cw[ii ][jj ] * p[ii - 1][jj ]) * cc[ii][jj];
        }
    }
    for (jj = 2; jj <= m_Nx2 + 1; jj += 2)
    {
        for (ii = 1; ii <= m_Nx1 + 1; ii += 2)
        {
            y[ii][jj] = (cs[ii ][jj + 1] * p[ii ][jj + 1]
                         + cw[ii + 1][jj ] * p[ii + 1][jj ]
                         + cw[ii ][jj ] * p[ii - 1][jj ]
                         + cs[ii ][jj ] * p[ii ][jj - 1]) * cc[ii][jj];
        }
    }
    // Do this multiplication for points (0,jj): west boundary of first level black points.
    ii = 0;
    for (jj = 1; jj <= m_Nx2 + 1; jj += 2)
    {
        y[ii][jj] = (cs[ii ][jj + 1] * p[ii ][jj + 1]
                     + cw[ii + 1][jj ] * p[ii + 1][jj ]
                     + cs[ii ][jj ] * p[ii ][jj - 1]) * cc[ii][jj];
    }
    // Do this multiplication for points (ii,0): south boundary of first level black points
    jj = 0;
    for (ii = 1; ii <= m_Nx1 + 1; ii += 2)
    {
        y[ii][jj] = (cs[ii ][jj + 1] * p[ii ][jj + 1]
                     + cw[ii + 1][jj ] * p[ii + 1][jj ]
                     + cw[ii ][jj ] * p[ii - 1][jj ]) * cc[ii][jj];
    }


    // Calcalate the matrix-vector multiplication y = D^-1 S p.

    // Do this multiplication for points (ii,jj): all first level red points outside boundary.
    for (jj = 1; jj <= m_Nx2 + 1; jj += 2)
    {
        for (ii = 1; ii <= m_Nx1 + 1; ii += 2)
        {
            q[ii][jj] = cc[ii ][jj ] * p[ii ][jj ]
                - cs[ii ][jj + 1] * y[ii ][jj + 1]
                - cw[ii + 1][jj ] * y[ii + 1][jj ]
                - cs[ii ][jj ] * y[ii ][jj - 1]
                - cw[ii ][jj ] * y[ii - 1][jj ];
        }
    }
    for (jj = 2; jj <= m_Nx2 + 1; jj += 2)
    {
        for (ii = 2; ii <= m_Nx1 + 1; ii += 2)
        {
            q[ii][jj] = cc[ii ][jj ] * p[ii ][jj ]
                - cs[ii ][jj + 1] * y[ii ][jj + 1]
                - cw[ii + 1][jj ] * y[ii + 1][jj ]
                - cs[ii ][jj ] * y[ii ][jj - 1]
                - cw[ii ][jj ] * y[ii - 1][jj ];
        }
    }
    // Do this multiplication for point (0,0): southwest corner of first level red points.
    ii = 0;
    jj = 0;
    q[ii][jj] = cc[ii ][jj ] * p[ii ][jj ]
        - cs[ii ][jj + 1] * y[ii ][jj + 1]
        - cw[ii + 1][jj ] * y[ii + 1][jj ];
    // Do this multiplication for points (ii,0): south boundary of first level red points.
    for (ii = 2; ii <= m_Nx1 + 1; ii += 2)
    {
        q[ii][jj] = cc[ii ][jj ] * p[ii ][jj ]
            - cs[ii ][jj + 1] * y[ii ][jj + 1]
            - cw[ii + 1][jj ] * y[ii + 1][jj ]
            - cw[ii ][jj ] * y[ii - 1][jj ];
    }
    // Do this multiplication for points (0,jj): west boundary of first level red points.
    ii = 0;
    for (jj = 2; jj <= m_Nx2 + 1; jj += 2)
    {
        q[ii][jj] = cc[ii ][jj ] * p[ii ][jj ]
            - cs[ii ][jj + 1] * y[ii ][jj + 1]
            - cw[ii + 1][jj ] * y[ii + 1][jj ]
            - cs[ii ][jj ] * y[ii ][jj - 1];
    }

} // end of matvec


